function Uint = unstr_int(x,grid,base,U,varargin)
% Uint = unstr_int(x,grid,base,U)
% Uint = unstr_int(x,grid,base,U,missing)
%
% Interpolation of the FE function U on the unstructured points x
% (passed as column vectors).
%
% First a search arlgorithm is used to find which element each point
% belongs to, and then the finite element basis is evaluated. The FE
% function U must be scalar and represented as a discontinuous field
% (i.e. as a two index array).

 np = size(x,2);

 xi = zeros(grid.d,np);
 Ue = zeros(base.pk,np);
 iei = 1;
 if(nargin>4)
   for i=1:np
     [iei,xii] = locate_point(x(:,i),grid,iei,1);
     if(iei>0)
       Ue(:,i) = U(:,iei);
       xi(:,i) = xii(2:grid.d+1);
     else
       Ue(:,i) = 0;
       xi(:,i) = 1/(grid.d+1);
       Ue(1,i) = varargin{1}/ev_pol(base.p_s{1},xi(:,i));
       iei = 1;
     end
   end
 else
   for i=1:np
     [iei,xii] = locate_point(x(:,i),grid,iei);
     Ue(:,i) = U(:,iei);
     xi(:,i) = xii(2:grid.d+1);
   end
 end

 PHI = zeros(base.pk,np);
 for i=1:base.pk
   PHI(i,:) = ev_pol(base.p_s{i},xi);
 end
 
 Uint = zeros(np,1);
 for i=1:np
   Uint(i) = PHI(:,i)'*Ue(:,i);
 end

return
